function ScalarFlux = mocsolver( track, BC, Q, SigmaT, SigmaS, MaxIt, EpsPhi )
% function ScalarFlux = mocsolver( track, BC, Q, SigmaT, SigmaS, MaxIt, EpsPhi )
%
% Inputs:
%       track           -- track information rom tracker
%       BC              -- boundary condition; 0=vacuum, 1=reflective (not done)
%       SigmaT          -- region-wise total cross-section
%       SigmaS          -- region-wise scattering cross-section
%       MaxIt           -- maximum inners
%       EpsPhi          -- relative tolerance for inner convergence
%
% Outputs:
%       ScalarFlux      -- the region-wise scalar flux
%
% References:  E.S.'s code inspired some of the structure
%
% J. Roberts, 03/05/2011
%

    % Preallocating average angular flux arrays (p,m,k,z)
    vol = zeros(2,1);
    % Preallocate the (previous) average angular flux.  The fourth
    % dimension if for indidual tracks.  Here, we preallocate space for the
    % maximum number of tracks.  This will yield a lot of zeros in the end,
    % and a more robust treatment would try to put the flux into fewer
    % dimensions using more complicated indexing (as in Hebert).  We stick
    % to this approach to ease the learning process.
    AngularFlux0 = zeros ( track.NumPolar, 2*track.NumAzimuth, ...
                           track.NumSpace, max(max(track.NumSeg)) ) ;
    % Preallocate current average angular flux .     
    AngularFlux  = AngularFlux0 ;
    
    % Note, the azimuth dimension is TWICE the number of azimuths tracked,
    % since we cover 0 to pi in the tracks, and the tracks from pi to 0
    % simply go in reverse along the first set.
    
    
    % Old and new (region-wise) scalar fluxes
    ScalarFlux0 = zeros( length(SigmaT), 1 ); 
    ScalarFlux  = ScalarFlux0;
    
    % Initialize the incident boundary fluxes and outgoing boundary fluxes
    IncidentAngularFlux = zeros(track.NumPolar,...
                                2*track.NumAzimuth,track.NumSpace); 
    ExitingAngularFlux  = IncidentAngularFlux; 
    
    % Initialize isotropic source
    S = Q/(4*pi);
    
    % Get the polar quantities
    [ SinTheta , WeightTheta ] = polar(track.NumPolar);
    
    Converged = false; % convergence flag  
    
    t1 = tic;
    
    % Running inner iterations
    for n = 1:1
             
        % Looping over polar divisions
        for p = 1:track.NumPolar

            % Looping over azimuths from 0 to Pi, then Pi to 2*Pi
            for mm = 1:track.NumAzimuth*2
                
                % mm is the FLUX index; m is the TRACK index
                if mm <= track.NumAzimuth
                    m = mm;
                else
                    m = mm-track.NumAzimuth;% 2*track.NumAzimuth - mm + 1; % e.g. mm = num+3 --> m = 2num-num-3 = num-3
                end

                if (m==0) 
                    a=1;
                end
%                 % Looping over tracks
%                 if mm <= track.NumAzimuth
%                     k1 = 1; k2 = track.NumSpace; dk = 1;
%                 else
%                     k1 = track.NumSpace; k2 = 1; dk = -1; % go along tracks in reverse
%                 end
                for k = 1:track.NumSpace

                    psi1 = IncidentAngularFlux(p,mm,k); % Angular flux at track entrance
                    
                    for z = 1:track.NumSeg(m,k) % loop over all segments
                        
                        if mm <= track.NumAzimuth
                            zz = z;
                        else
                            zz = track.NumSeg(m,k)-z+1;
                        end
                        
                        % get my material
                        if (mm==11)
                            lala=1;
                        end
                        i = track.SegData(m,k).SegRegIdx(zz);
                        
                        % Updating flux at segment entrance
                        psi0 = psi1 ;
                        
                        % Checking for transparency (to avoid division by zero )
                        if SigmaT(i) == 0
                            psi1 = psi0 ;
                            AngularFlux(p,mm,k,zz) = psi1 ; % Not sure if this is right
                        else                        
                            % Many of these values (A,B,C,etc) should be
                            % able to be precomputed.  
                            t    = track.SegData(m,k).SegLength(zz)/SinTheta(p);    % Polar-corrected segment length
                            A    = exp(-SigmaT(i)*t);                              % A, B, and C
                            B    = (1.0-A)/SigmaT(i);                              % from Hebert's book
                            C    = t/SigmaT(i) * (1 - (1-A)/(t*SigmaT(i)));
                            psi1 = A*psi0 + B*S(i);                % Angular flux at segment exit
                            if (~(psi1 >= 0))
                                lala==1;
                            end
                            AngularFlux(p,mm,k,zz) = ...
                                (1/t)*(B*psi0+C*S(i));             % Average flux along this segment <== Did I correct Vujic's typo??
                            %AngularFlux(p,mm,k,zz)=1; psi1=0;
                            
                            vol(i) = vol(i) + AngularFlux(p,mm,k,zz) ;%* track.Weight(m) *pi;
                        end
                        
                    end % segments
                    
                    ExitingAngularFlux(p,mm,k) = psi1; % exiting flux is the last one computed
                    
                end % space
                
            end % azimuth
            
        end % polar
                
        disp('summmed psi'), vol
        % Update scalar flux and check for convergence
        ScalarFlux = scalar(track,AngularFlux); 
        FluxErr    = (ScalarFlux-ScalarFlux0)./ScalarFlux ;  % Scalar flux fractional error vector
        if norm ( abs ( FluxErr ), inf ) < EpsPhi % Converged ( using infinity norm )
            Converged = true;
            break
        else % Go to next iteration
            % Update Source
            S = (Q + ScalarFlux.*SigmaS)/(4.0*pi);
            % Update Boundary Conditions
            %sbound = calcsbound (track,ExitingAngularFlux,BCtype);            
            % Store most recent scalar flux
            ScalarFlux0 = ScalarFlux;
        end
        
       % Reporting iteration completion
        fprintf ('Completed inner iteration %g\n',n);
        fprintf ('Cumulative time %f\n',toc(t1));
        
    end % for over iterations
    
    if (Converged)
        fprintf ('Inner iterations converged successfully in %g iterations .\n',n);
    else
        disp('Failure to converge inner iterations !');
    end    

end % function mocsolver


function ScalarFlux = scalar(track,AngularFlux)
% function scalar(track,AngularFlux)
% 
% This function computes region-wise scalar fluxes from the tracking
% information and the full angular flux. 
%
%   Inputs:
%           track       -- tracking informatin
%           AngularFlux -- angular flux in angle, space, and segment, i.e.
%                          [NumPolar,NumAzimuthal,NumSpace,MaxNumSegments]
%   Output:
%           ScalarFlux  -- region scalar fluxes
%
    % track - Tracking information data structure
    % flux - Average angular flux array (p,m,k,z)
    % sflux - Scalar flux vector (i)
    
    % Get the polar quantities
    [ SinTheta, WeightTheta ] = polar(track.NumPolar);
   
    % Preallocating arrays
    ScalarFlux          = zeros(track.NumReg,1);                 % Scalar flux vector (i)
    PolarIntAngularFlux = zeros(track.NumAzimuth*2,track.NumReg);  % Polar integrated, region-average angular flux
    Volume              = zeros(track.NumAzimuth*2,track.NumReg);  % Polar integrated, region-average angular flux
    
    % Computing average angular fluxes ( averaging over volumes ), Psi_{im}
    % in eq. 9 o Vujic.  Defined by
    %   Psi_{im} = Sum( Psi_{imk}*s_{imk}*del_{k} ) / 
    %              Sum( s_{imk}*del_{k} )
    % but the del_{k} is constant and cancels.  
    for p = 1:track.NumPolar
        
        wp = 2*WeightTheta(p); % double since int(dmu,-1,1)=2
        
        for mm = 1:track.NumAzimuth*2
            
            % mm is the FLUX index; m is the TRACK index
            if mm <= track.NumAzimuth
                m = mm;
            else
                m = mm-track.NumAzimuth;%2*track.NumAzimuth - mm + 1;
            end
            
            % Looping over tracks
%             if mm <= track.NumAzimuth
%                 k1 = 1; k2 = track.NumSpace; dk = 1;
%             else
%                 k1 = track.NumSpace; k2 = 1; dk = -1; % go along tracks in reverse
%             end
            for k = 1:track.NumSpace
                
                for z = 1:track.NumSeg(m,k)
                    
                    if mm <= track.NumAzimuth
                        zz = z;
                    else
                        zz = track.NumSeg(m,k)-z+1;
                    end
                    
                    % get my material
                    i = track.SegData(m,k).SegRegIdx(zz);
                     
                    % add track contribution to the region
                    PolarIntAngularFlux(mm,i) = PolarIntAngularFlux(mm,i) ...
                        + wp * AngularFlux(p,mm,k,zz) * ...
                        track.SegData(m,k).SegLength(zz);%/SinTheta(p);
                    
                    if ( ~(PolarIntAngularFlux(mm,i)>=0) && ~(PolarIntAngularFlux(mm,i)<=0) )
                        lala=1;
                    end
                    % add contribution to region volume
                    Volume(mm,i) = Volume(mm,i) ...
                        + wp * track.SegData(m,k).SegLength(zz);%/SinTheta(p);                        
                        
                end % segment
            
            end % space
            
        end % azimuth
        
    end % polar
    %PolarIntAngularFlux = PolarIntAngularFlux ./ Volume;
    
    % Integrating over azimuthal angle
    for i = 1:track.NumReg
        sf = 0;
        for mm = 1:track.NumAzimuth*2
            if mm <= track.NumAzimuth
                m = mm;
            else
                m = mm-track.NumAzimuth;%2*track.NumAzimuth - mm + 1;
            end
            sf = sf + 1.0*pi*track.Weight(m)*PolarIntAngularFlux(mm,i)*track.PertDist( m ) / track.Volume(i);
        end
        ScalarFlux(i) = sf;
    end
    
end